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Executive  summary 


Objectives 

To  assess  the  feasibility  of  bridging  the  gap  between  atomistic  scale  and  microstructural 
scale  modelling. 

Why 

If  we  can  understand  how  variations  on  the  atomistic  level  modify  the  microstructure,  we 
will  be  in  a  position  to  engineer  microstructures  by  making  modifications  at  the  atomic 
level. 

Approach 

The  possible  approaches  to  the  generation  of  parsimonious  rules  for  the  Cellular 
Automata  (CA)  engine  from  atomistic  scale  simulations  were  first  reviewed  (see 
preliminary  report)  and  here  we  demonstrate  the  feasibility  by  addressing  a  specific 
problem. 

Conclusions 

1.  CA  can  be  used  to  simulate  microstructural  evolution  using  rules  generated  by 
atomistic  simulations. 

2.  This  was  demonstrated  using  the  specific  example  of  the  adsorption/desorption  of 
inert  gas  atoms  on  a  Ca  (1 1 1)  surface. 

3.  The  computational  benefits  if  using  atomistic  based  rules  cellular  automata  (AR- 
CA)  are  considerable.  We  have  shown  a  computational  gain  over  MD  simulations 
of  the  same  size  of  10^. 

4.  The  computational  gain  can  be  increased  significantly  by  using  multiple  atom  CA 
cells  of  a  hierarchical  approach. 

5.  Specific  conclusions  from  the  first  implementation  of  an  AR-CA  model  are: 

i.  island  growth  is  observed; 

ii.  island  growth  is  a  strong  fimction  of  temperature  and  gas  atom  species; 
and 

iii.  island  growth  rates  depend  on  the  radius  of  the  boundaries  -  consequently, 
a  ripening  process  of  island  growth  is  observed  over  a  specific  temperature 
range. 
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Simulating  the  Evolution  of  an  Atomic  Layer  of 
Neon  or  Argon  on  a  Calcium  (111)  Surface  Using 
Atomistic  Based  Cellular  Automata  Rules 

1  Motivation 

Cellular  automata  (CA)  methods  are  computationally  simple  and  as  such  are  attractive 
alternatives  to  Molecular  Dynamics  (MD)  and  Monte  Carlo  (MC)  methods  for  simulating 
the  dynamic  evolution  of  solid  state  systems.  This  is  particularly  true  in  the  present 
context  where  we  wish  to  simulate  microstructural  evolution  at  the  atomic  level.  As  will 
become  clear,  the  size  of  the  growth  structures  associated  with  the  gas  atom  layers 
requires  us  to  consider  40,000  gas  atom  sites.  Furthermore,  we  find  that  microstructure 
can  evolve  slowly  with  respect  to  the  association  and  evaporation  of  gas  atoms  which 
implies  that  many  tens  of  thousands  of  atomic  vibrations  would  have  to  be  simulated 
explicitly  if  an  MD  method  were  to  be  employed.  This  would  be  a  major  undertaking  for 
an  MD  code  whereas  using  CA  it  can  be  carried  out  in  a  matter  of  minutes  on  a  machine 
of  modest  power. 

Of  course,  the  system  we  vdll  consider  here  is  rather  specific  and  not  itself  of  great 
technological  importance.  However,  to  the  authors  knowledge,  this  is  the  first  time  that 
atomistic  simulation  data  has  been  used  as  the  basis  for  CA  rules.  As  such,  this  study 
presents  a  powerful  new  methodology  for  bridging  between  atomistic  and  microstructural 
scale  simulations. 

2  The  Problem 

We  will  describe  the  evolution  of  a  single  layer  of  neon  or  argon  on  a  calcium  (111) 
surface.  This  is  treated  as  a  two  dimensional  problem,  thus  it  is  assumed  that  if  gas  atoms 
do  form  a  second  covering  layer,  this  is  not  of  any  significance  to  the  evolution  of  the  first 
layer. 

Since  this  is  an  atomic  scale  problem,  it  is  necessary  to  account  explicitly  for  the 
interaction  between  a  gas  atom  and  the  calcium  surface  and  between  a  gas  atom  and  its 
neighbouring  gas  atoms  (inter  gas  atom  interactions).  Interatomic  potentials  were 
therefore  derived  which  describe  the  interaction  energy  between  atomic  species  as  a 
function  of  interatomic  separation  (these  will  be  described  in  detail  below).  Thus,  once 
we  know  the  position  of  a  gas  atom  on  the  surface  and  its  number  of  neighbouring  gas 
atoms,  it  is  straightforward  to  determine  its  stability.  In  the  approach  used  here,  these 
energies  are  translated  into  local  rules  for  the  CA  engine. 

The  starting  point  for  all  the  simulations  described  here  is  a  bare  metal  surface.  The  gas 
atom  layer  is  deposited  onto  the  surface,  from  the  gas  phase,  as  a  flux  whose  magnitude 
can  be  varied.  Gas  atoms  may  also  desorbe  from  the  surface.  The  surface  therefore 
evolves  via  a  combination  of  adsorption  and  desorption  processes  both  of  which  are 
temperature  dependent  via  the  CA  rules. 
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Clearly  it  is  essential  to  understand  the  physical  basis  of  the  CA  rules.  Thus  in  the  next 
section  we  will  consider  the  (111)  surface  structure  and  the  positions  of  the  gas  atoms  in 
detail. 


2.1  Crystallography 

Calcium  metal  has  an  FCC  structure  so  that  the  (1 1 1)  surface  exhibits  a  close  packing  of 
Calcium  atoms.  These  may  be  regarded  as  "A"  sites.  When  a  gas  atom  attaches  itself  to 
the  surface  it  will  sit  m  a  site  which  is  formed  by  three  metal  atoms  since  this  maximizes 
its  interaction  with  the  metal  surface.  It  is  clear  from  the  figure  below  that  two  sets  of 
interstitial  sites  exist,  these  can  be  regarded  as  "B"  sites  and  "C"  sites. 

gas  atom 
at  a  "B"  site 


If  we  examine  the  (111)  surface  more  closely  it  is  clear  that  the  distance  between  an 
interstitial  "B"  site  and  a  nearest  neighbouring  "C"  site  is  rather  short.  Indeed  the  gas 
atom  in  the  figure  above  (which  represents  neon  in  scale  to  calcium)  overlaps 
considerably  onto  the  "C"  site.  As  such,  if  a  gas  atom  resides  in  a  "B"  site,  it  is  not 
possible  for  a  second  gas  atom  to  occupy  an  nearest  neighbour  "C"  site  due  to  steric 
hinderence.  It  is  important  to  realise  that  this  is  true  only  because  of  the  specific  matching 
of  the  size  of  the  Ca  (1 1 1)  surface  to  the  sizes  of  the  Ar  and  Ne  ions. 

Let  us  now  consider  the  evolution  of  the  gas  atom  layer  as  follows.  When  the  first  few 
gas  atoms  impinge  on  the  surface,  they  may  reside  in  either  A  or  B  sites.  As  more  atoms 
collect,  patches  or  islands  of  B-type  and  C-type  atoms  form.  In  the  diagram  below,  it  is 
clear  that  the  distance  between  two  nearest  neighbour  B-type  gas  atoms  is  the  same  as 
between  two  C-type  gas  atoms.  However,  the  distance  between  allowed,  none  hindered 
B-type  and  C-type  atoms  is  greater.  This  results  in  an  area  of  surface  with  a  lower  density 
of  gas  atoms  which  will  be  termed  a  growth  fault.  In  three  dimensions,  it  would  give  rise 
to  a  stacking  fault.  This  is  important  to  the  stability  of  a  gas  atom  on  the  surface  since 
when  gas  atoms  interact  over  a  longer  distance,  their  interaction  energy  is  lower.  Indeed, 
this  is  the  driving  force  which  results  in  the  formation  and  growth  of  islands  of  gas  atoms 
which  are  occupying  only  one  type  of  site. 
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Growth 

Fault 


Growth 

Fault 


2.2  Interatomic  potential  functions 

Inert  gas  atoms  interact  through  two  terms,  electron  cloud  overlap  repulsion  and  van  der 
Waals  attraction.  In  this  study,  these  are  manifest  in  Lennard-Jones  potentials  where,  if  pj 
is  the  distance  between  atoms  i  and  j,  the  interaction  energy,  Ejj,  is  given  by: 

Eij  =  Arij-'"-Br,-‘ 

where  A  and  B  are  constants  specific  to  atoms  i  and  j.  For  each  potential  type,  the 
constant  "B"  corresponds  to  the  van  der  Waals  interaction  and  can  be  determined  from  the 
polarisabilities  of  the  interacting  ions  via  the  Slater  Kirkwood  formulae  (see  Fowler  et  al. 
Mol.  Phys.,  56,  83,  (1985)).  For  the  inert  gases,  the  "A"  constants  were  determined  by 
ensuring  that  the  appropriate  potential  reproduced  the  lattice  parameter  of  the  inert  gas 
lattice  (note  that  both  neon  and  argon  exhibit  FCC  latticies).  For  the  inert  gas  -  metal 
interaction,  the  "A"  parameter  was  chosen  so  that  the  potential  has  a  minimum 
corresponding  to  ry  =  gas  radius  +  metal  radius. 

Table  I:  Interatomic  Potential  Parameters 


Potential  type 

A  (eVAng^^) 

B  (eVAng*) 

Ca-Ne 

12276.6 

12.54 

Ca-Ar 

64447.4 

40.90 

Ne-Ne 

1953.89 

4.14 

Ar  -  Ar 

53845.3 

40.90 

Once  the  parameters  have  been  determined,  it  is  possible  to  calculate  the  interaction 
energies,  which  define  the  stability  of  an  atom  at  a  specific  site  on  the  surface.  Usually 
with  atomistic  simulation,  a  calculation  would  be  carried  out  during  which  all  atoms 
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would  be  allowed  to  relax  to  zero  force.  In  ionic  systems,  where  long  range  Coulomb 
forces  operate,  this  is  critical.  However,  in  this  highly  symmetric  system,  where  energies 
are  described  by  short  range  van  der  Waals  interactions,  a  more  simple  approach  is 
justifiable.  First,  the  interaction  of  a  gas  atom  with  the  substrate  is  given  by  three  times 
the  energy  of  interaction  of  a  gas  atom  with  a  single  calcium  atom  at;  nj  =  gas  radius  + 
metal  radius.  Second,  the  interaction  energies  between  nearest  neighbour  "B"  gas  atoms  is 
given  by  the  value  of  the  potential  at;  ry  =  Ca  -  Ca  nearest  neighbour  distance  =  twice  the 
Ca  metal  ion  radius.  Remember,  the  separation  of  the  gas  atoms  is  defined  by  the 
substrate.  The  same  is  true  for  interactions  between  nearest  neighbour  "C"  gas  atoms.  The 
total  interaction  between  a  gas  atom  and  its  nearest  neighbours  depends  on  the  number  of 
nearest  neighbours.  This  is,  of  course,  a  key  feature  of  the  model. 

A  further  simplification  is  made  to  describe  the  interactions  between  "B"  site  gas  atoms 
and  nearest  allowed  "C"  site  gas  atoms:  they  are  assumed  to  be  zero.  Although  this  may 
seem  a  somewhat  draconian  assumption,  the  increase  in  the  gas  -  gas  distance  in  the  (2,- 
1,-1,0)  direction  (hexagonal  coordinate  system)  over  the  nearest  "B"  -  "B"  is  1.5366 
which  implies  a  reduction  in  the  gas  -  gas  interaction  by  a  factor  of  13.  The  reduction  in 
interaction  in  the  (0,1, -1,0)  direction  is  somewhat  less,  a  factor  of  2.2.  Given  this,  it 
would  be  interesting,  in  future  work,  to  account  for  the  "B"  -  "C"  interaction  in  the  (0,1,- 
1,0)  direction  which  may  lead  to  some  anisotropic  growth  effects. 

Given  the  above  rules,  it  is  simple  to  calculate  the  energy  that  a  gas  atom  would 
experience  if  it  were  to  occupy  a  specific  surface  site.  This  will  be  known  as  the 
attachment  energy.  Eg.  Note  that  because  of  the  way  that  the  interaction  potentials  are 
defined,  negative  attachment  energies  are  favourable. 

This  leaves  us  with  the  problem  of  relating  the  attachment  energy  to; 

i) .  The  probability  that  a  fi-ee  gas  atom,  incident  on  the  surface,  will  become  attached  to 
the  surface  and 

ii) .  That  over  a  given  time  period,  an  attached  gas  atom  will  desorbe  fi-om  the  surface. 

This  will  be  the  subject  of  the  next  two  sections. 


2.3  Adsorption 

The  probability  for  adsorption  is  the  fraction  of  incoming  atoms  that  become  attached  to 
the  surface,  ?att-  The  incoming  atoms  are  assumed  to  have  a  Maxwell  distribution  of 
kinetic  energies,  e; 

n(e)  =  (l/kT)exp(-e/kT) 

The  jfraction  of  incoming  atoms  that  stick  to  the  surface  depends  on  the  substrate's  ability 
to  dissipate  a  sufficient  proportion  of  the  kinetic  energy  of  a  colliding  atom  for  the  atom 
to  become  trapped  in  the  potential  well.  Thus,  we  make  the  assumption  that  any  atom 
with  kinetic  energy  less  than  Ea  will  stick.  Thus; 


Page  4 


Final  Report  BOARD  contract  F61708-96-W0300 


Sept.  1997 


^ati  ~ 


Ed  j 

J  j^exp(-e/kT>3fe 


=  l-exp(-Ea/kT) 


2.4  Desorption 

Assuming  a  single  step  activated  process,  the  probability,  that  desorption  occurs 
between  t  and  i  +  dt\s  given  by; 


Prfe^CRd,!)  =  Rdexp(-Rdt) 
where  the  desorption  rate,  Rd,  is  given  by; 

Rd  =  gdf;,exp(Ea/kT) 

where  T  is  temperature  in  Kelvin,  k  is  Boltzmann's  constant,  gd  is  a  geometric  factor 
(which  we  take  to  be  unity  in  this  case)  and  fo  is  the  attmept  frequency.  Thus,  the  total 
number  of  atoms  which  desorbe  over  a  single  time  step,  st  is; 

Prfe^  =  jRdexp(-Rdt)it  =  l-exp(-RdSt) 

0 

The  attempt  frequency  is,  to  a  good  approximation,  equivalent  to  the  vibrational 
frequency  of  the  gas  atom  on  the  surface.  This  can  be  determined  from  the  interatomic 
potentials.  Here  the  energy  of  each  gas  atom  was  determined  as  a  function  of  its  distance 
perpendicular  to  the  surface,  assuming  a  frill  compliment  of  neighbouring  gas  atoms  (i.e. 
a  symmetric  site).  Harmonic  functions  were  then  fitted  to  each  energy  profile  to  yield  the 
force  constants.  The  characteristic  frequencies  were  then  determined  from  the  force 
constants;  these  were  4.45x10^°  s'*  for  Ar  and  4.59x10'*  s'*  for  Ne. 

If  the  attempt  frequency  is  not  known,  it  can  be  assumed  to  be  equal  to  the  reciprocal  of 
the  time  step,  in  which  case,  the  desorption  probability  simplifies  to; 

Pdfei  =  1  -  exp{(exp-Ea/kT)} 

which  is  approximately  equal  to  exp(-Ea/kT)  when  the  value  of  Ea/kT  is  small.  However, 
if  this  approximation  is  made,  we  caimot  relate  the  flux  to  any  real  time  scale  and  the  only 
simulations  that  are  definitive  are  those  that  run  to  equilibrium. 


2.5  The  flux 

In  this  work,  the  flux  is  defined  to  be; 
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In  - 


(number  of  atoms  striking  the  surface) 
(number  of  surface  sites)(time  step) 


This  is  an  atomistic  model  which  naturally  lends  itself  to  the  implementation  of  the  type 
of  relationship  discussed  above.  The  flux  used  here  was  15%,  that  is,  the  number  of  gas 
atoms  which  were  allowed  to  interact  with  the  surface  corresponded  to  15%  of  the  total 
number  of  surface  sites  (B  and  C).  It  is  important  to  realise  that  some  of  these  atoms 
could  attempt  to  interact  with  the  same  surface  site.  In  this  case,  only  one  atom  is  allowed 
to  be  successful.  This  way  of  implementing  the  flux  is  realistic  since  gas  atoms  interact 
with  the  surface  continuously  rather  than  in  the  step  wise  discrete  manner  of  CA. 
However,  if  two  atoms  attempt  to  simultaneously  occupy  nearest  neighbour  sites,  none 
are  allowed  to  be  successful.  This  is  not  physically  realistic  for  the  long  CA  time  steps 
used,  but  it  is  computationally  much  simpler  than  any  alternative  choice  algorithm  yet 
discovered. 


Once  the  flux  value  is  chosen,  for  a  given  box  size,  this  defines  the  number  of  gas  atoms, 
each  of  which  are  assigned  a  random  number  from  1  to  the  total  number  of  surface  sites. 
Integer  division  by  2  is  then  carried  out  to  yield  an  integer  and  a  remainder.  If  the 
remainder  is  zero,  the  atoms  will  strike  a  "B"  site,  if  it  is  one,  it  will  strike  a  "C"  site.  The 
integer  value  assigns  the  atom  to  a  surface  cell  which  contains  a  "B"  site  and  a  "C"  site. 
The  number  of  surface  cells  is,  of  course,  half  the  number  of  surface  sites. 

It  would  be  a  great  advantage  to  be  able  to  relate  the  flux  to  the  pressure  of  the  gas  above 
the  surface.  In  this  case,  we  can  define  the  flux  in  terms  of  the  number  of  atoms  striking  a 
surface  area  in  unit  time  and  relate  this  to  the  velocity  of  an  ideal  gas. 

(number  of  atoms  striking  the  surface)  pc  "■ 

(surface  area)(time  step)  kT 

where  p  is  the  partial  pressure  of  the  gas  and  c*  is  the  mean  velocity  of  the  incoming 
atoms.  Then, 


so  that  the  flux  can  be  related  to  the  pressure  by 

2d,UMkT7 

where  ds  is  the  surface  density  of  absorption  sites  and  M  is  the  mass  of  the  absorbing 
species. 
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3  The  CA  rule  set 

To  a  great  extent,  the  CA  rule  set  is  defined  by  the  ideas  introduced  in  the  preceding 
sections.  These  can  now  be  presented  formally. 

3.1  Site  selection  or  flux  rules 

These  concern  the  assignment  to  the  incident  atoms,  of  a  single  random  number  which 
determines  which  surface  cell  the  atom  will  strike  and  weather  that  will  be  to  the  "B"  or 
"C"  site  within  that  cell. 

The  first  rule  assigns  to  each  atom  incident  on  to  the  surface  a  random  number  from  one 
to  the  total  number  of  surface  sites. 

This  number  then  undergoes  integer  division  by  two. 

The  second  rule  uses  the  value  of  the  remainder  fi-om  the  integer  division.  If  the 
remainder  is  zero,  the  atom  is  incident  on  to  a  "B"  site,  if  the  remainder  is  one,  the  atom 
is  incident  on  to  a  "C"  site. 

The  third  rule  uses  the  integer  value  fi-om  the  integer  division.  The  integer  value  assigns 
the  atom  to  a  particular  cell  on  the  surface. 

The  combination  of  a  surface  cell  and  the  assignment  of  "B"  or  "C"  site  yields  a  unique 
interstitial  site. 

3.2  Unavailable  site  rules 

The  first  rule  concerns  occupied  sites.  When  the  flux  directs  a  gas  atom  to  an  interstitial 
site  that  is  already  occupied,  the  gas  atom  will  not  become  attached  to  the  surface. 

The  second  rule  concerns  the  hinderence  that  an  occupied  "C"  site  has  to  "B"  site 
occupation  and  vice  versa.  When  the  flux  directs  a  gas  atom  to  a  "B"  interstitial  site  that 
has  a  gas  atom  at  a  nearest  neighbour  "C"  site,  the  gas  atom  will  not  become  attached  to 
the  surface.  There  is  an  equivalent  hinderence  rule  for  gas  atoms  directed  to  a  "C"  site. 

3.3  Available  site  attachment  rule 

This  rule  concerns  the  probability  of  a  gas  atom  becoming  attached  once  it  has  been 
determined  that  a  site  is  available  for  occupation.  It  assumes  that  the  attachment  energy, 
Ea  has  been  calculated.  The  probability,  Pa«(Ea),  that  a  gas  atom  will  become  attached  to 
the  site,  to  which  it  has  been  directed  to  by  the  flux,  is  given  by  a  Bolzmann  factor; 


Pa«(Ea)  =  1  -  exp(Ea/kT) 
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3.4  Detachment  rule 

Within  the  approximations  discussed  above,  at  each  time  step,  all  attached  gas  atoms  are 
given  the  opportunity  to  desorbe.  The  probability  that  a  gas  atom  will  desorbe  from  the 
surface,  P(fei(Rd)  is  given  by  a  Bokmann  factor; 

P*^(Rd,St)  =  1  -  exp(-R<iSt) 


4  Visualization 

To  perform  the  calculations,  the  program  uses  a  diamond  portion  of  the  hexagonal  lattice. 
In  the  results  discussed  here,  this  was  chosen  to  be  200  by  200  sites,  i.e.  40,000  sites  in 
total.  Each  pair  of  sites  on  the  surface  is  indexed  by  coordinates  i  and  j  as  shown  in  the 
figure  below.  Each  site  in  a  pair  is  indexed  by  a  1  or  2.  In  order  to  minimize  memory 
usage,  the  results  of  the  calculations  are  shown  on  a  square  grid  rather  than  a  diamond 
grid,  allowing  one  pixel  to  represent  one  cell.  The  colour  red  indicates  that  site  1  of  cell 
(i,j)  is  filled,  the  collar  blue  indicates  that  site  2  is  filled,  and  the  collar  black  indicates  that 
neither  site  is  filled. 
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5  Results 

The  following  are  typical  structures  evolved  for  a  series  of  temperatures  and  fluxes  for 
both  argon  (left  hand  side)  and  neon  (right  hand  side)  on  a  calcium  (111)  surface.  These 
results  may  be  viewed  as  mpeg  movies  in  the  electronic  version 
(http://abulafia.mt.ic.ac.uk/USAF/). 
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The  first  implication  of  these  simulations  is  that  for  Ne  and  Ar,  the  rate  of  island  growth 
is  highly  temperature  dependent.  Second,  the  temperature  at  which  growth  is  a  maximum 
is  ~35K  or  Ne  and  ~130K  for  Ar.  The  increase  in  growth  rate  with  temperature  is  clearly 
due  to  the  increase  in  the  kinetics.  Presumably  if  the  lower  temperature  simulations  were 
implemented  over  more  time  steps,  more  pronounced  island  growth  would  eventually 
occur.  However,  at  very  low  temperatures,  this  would  take  an  inordinately  long  time. 
Finally,  the  limited  island  growth  beyond  the  temperature  of  maximum  growth  rate  is  due 
to  the  instability  of  islands  due  to  rapid  evaporation. 

In  the  next  section  we  shall  use  these  results  to  formulate  an  atomistic  explanation  for 
island  growth.  This  is  only  possible  because  the  CA  rules  are  atomistically  based.  The 
potential  importance  of  this  work  is  now  apparent.  If  we  can  understand  how  variations 
on  the  atomistic  level  modify  the  microstructure,  we  will  be  in  a  position  to  engineer 
microstructures  by  making  modifications  at  the  atomic  level. 
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5.1  Interpretation 

Consider  the  specific  interface  between  two  islands  depicted  below.  Above  each  atom  is 
its  calculated  attachment  energy.  It  is  clear  that  interface  sites  have  lower  attachment 
energies  compared  to  sites  which  have  a  full  complement  of  neighbours.  This  has  two 
implications.  First,  interface  atoms  have  a  higher  probability  of  evaporating.  Second, 
atoms  that  are  incident  on  a  interface  site  are  less  likely  to  become  attached. 
Consequently,  the  evolution  of  the  atomic  distribution  is  more  active  at  interfaces. 


Growth 

Fault 


Growth 

Fault 


However,  the  situation  is  complicated  by  the  fact  that  some  interface  sites  have 
substantially  lower  attachment  energies  than  others.  Thus  the  more  highly  convoluted  an 
interface,  the  more  quickly  it  becomes  modified  through  the  evaporation  attachment 
process.  Furthermore,  not  only  is  there  an  energetic  imperative  for  this  to  occur,  there  can 
also  be  a  crystallographic  hinderence  effect.  This  can  most  easily  be  understood  if  we 
consider  a  straight  (1,0, -1,0)  interface. 
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For  the  interface  to  progress  upwards,  expanding  the  "C"  island,  a  "B"  atom  must  be  lost. 
This  is  shown  in  the  next  figure. 


Growth 

Fault 


However,  both  of  the  potential  "A"  sites  are  hindered  by  "B"  type  atoms.  Thus  a  second 
"B"  atom  must  evaporate,  adjacent  to  the  first  site.  This  situation  is  shown  in  the  next 
figure. 


Growth 

Fault 


It  is  now  possible  for  a  single  "C"  type  atom  to  attach  itself  to  the  surface,  effectively 
expanding  the  "C"  island  as  shown.  Of  course,  the  probability  that  two  atoms  are  lost 
adjacent  to  each  other  is  correspondingly  small.  Furthermore,  this  "C"  type  atom  has  only 
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two  neighbouring  "C"  type  atoms  which  means  that  its  attachment  energy  is  lower  than  if 
a  "B"  type  atom  were  to  reattach  itself.  Thus  the  expansion  of  a  straight  boundary  is 
subject  to  energetic  and  probabalistic  barriers. 


Growth 

Fault 


By  running  the  argument  the  opposite  way  round  it  is  possible  to  see  that  these  type  of 
convoluted  structures  will  tend  to  be  removed  as  the  simulation  progresses.  Indeed,  by 
applying  this  analysis  to  other  situations,  one  can  conclude  that  boundaries  should 
straighten  on  an  atomistic  level.  Of  course,  this  is  a  well  known  experimental  observation 
which  is  usually  discussed  in  terms  of  surface  tension. 

At  temperatures  close  to  the  maximum  for  island  growth,  the  lower  attachment  energies 
of  interface  ions  can  have  the  effect  of  forming  regions,  adjacent  to  interfaces  which  are 
denuded  of  ions.  These  show  up  on  the  simulations  as  dark  areas.  However  these  are 
usually  small  and  hard  to  see. 

6  Possible  limitations  the  CA  model 

One  of  the  potential  problems  that  has  become  apparent  during  the  progress  of  this  study 
concerns  the  reliability  of  predicted  rates  of  island  growth.  Here  we  calculate  an  energetic 
advantage  of  atoms  in  the  interior  of  islands  over  atoms  at  an  interface  region  of  between 
10-50  %.  The  simulation  is  then  run  for  2,000  time  steps  which  corresponds  to 
approximately  a  microsecond.  If  the  relative  attachment  energies  were  slightly  different, 
the  predicted  temperatures  at  which  the  maximum  in  island  growth  occurs  could  be  quite 
different.  This  is  a  general  problem  which  should  be  given  more  attention.  However,  for 
many  systems,  experimental  data  may  be  available  in  which  case,  the  energies  could  be 
learnt  by  running  the  simulation  and  matching  this  simple  parameter. 

In  this  model,  each  CA  cell  contains  two  sites  of  which  only  one  may  be  occupied.  This  is 
fixed.  Although  we  are  able  to  simulate  a  greater  number  of  atoms  (~10^),  over  a  far 
longer  total  time  interval  (~10‘®  s)  than  any  previous  MD  type  simulation,  for  the 
modelling  of  a  three  dimensional  atomic  microstructure,  we  are  still  very  limited.  For 
example,  a  typical  grain  in  a  conventional  microstructure  would  be  at  least  ~0.1 
micrometres  across.  This  would  encompass  ~5xl0^  atoms.  Clearly  we  would  wish  to 
include  more  than  one  grain  in  any  simulation!  In  addition,  a  typical  heat  treatment  would 
last  for  at  least  10^  seconds.  The  approach  to  predicting  real  three  dimensional 
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microstructural  evolution  will  have  to  be  fundamentally  different.  The  most  obvious 
modification  would  be  to  include  multiple  atomic  sites  in  a  single  CA  cell.  Nevertheless, 
the  properties  of  atoms  in  the  CA  cell  could  and  probably  would  have  to  be  learnt  from  a 
combination  of  smaller  scale  CA  or  even  MD  simulations.  As  such,  we  envisage  a 
hierarchy  of  models  which  culminate  in  a  full  blown  real  scale  simulation.  In  this  regard, 
the  scales  of  the  present  simulations  are  important. 


7  Ideas  for  future  work 

7.1  Specific  to  the  present  code 

In  the  short  term,  it  would  be  very  useful  to  run  energy  minimisation  atomistic 
simulations  of  atoms  in  a  variety  of  situations  and  formulate  a  more  complete  look-up 
table  of  attachment  energies  thereby  formally  completing  the  atomistic  to  CA  loop. 

Certainly  a  number  of  tools  would  be  useful  in  providing  analysis  of  the  simulations.  The 
first  would  be  a  way  of  reporting  island  size.  This  would  allow  us  to  quantify  the  extent  of 
the  grain  growth  over  a  given  simulation  time  and  grain  growth  rates  as  a  function  of 
temperature  and  pressure. 

There  are  a  number  of  additional  topics  which  the  present  code  should  be  used  to  address. 
The  first  is  to  introduce  another  atom  type.  A  particularly  interesting  example  would  be 
xenon  which  due  to  its  greater  mass  and  higher  polarizability  would  stick  to  the  surface 
and  attract  neighbouring  Ar  or  Ne  atoms  strongly.  It  is  envisaged  that  xenon  would  then 
act  as  an  island  growth  inhibitor,  effectively  pinning  convoluted  structures.  Again  this  is 
not  a  technologically  important  example  but  it  would  expand  our  assessment  of  what  this 
type  of  model  can  simulate  successfully. 

In  these  simulations,  we  start  with  an  empty  surface.  As  an  alternative,  we  could  begin 
with  an  idealised  boundary  and  investigate  how  this  evolves.  For  example,  to  what  extent 
will  a  perfect  straight  (1,0,- 1,0)  boundary  become  convoluted  over  a  period  of  time  at  a 
given  temperature?  Certainly  it  should,  due  to  entropic  effects,  but  to  what  extent  is  less 
clear. 

It  would  also  be  possible  to  investigate  the  specific  rate  of  shrinkage  of  a  boundary  as  a 
function  of  island  radius  by  initiating  a  simulation  with  a  single  circular  island  of 
specified  radius. 

7.2  Extensions  to  the  present  atomistic  based  code 

It  would  be  a  simple  matter  to  extend  the  present  code  to  include  diffusion  since  in  this 
2D  context,  this  corresponds  to  activated  hops  into  vacant  sites.  The  activation  energy  can 
easily  be  calculated  using  energy  minimisation  atomistic  simulation.  Within  islands,  if 
migration  is  mediated  by  a  single  vacant  site,  this  corresponds  to  "B"  to  "B"  or  "C"  to  "C" 
jumps.  However,  at  a  boundary  between  "B"  and  "C"  islands,  there  is  the  possibility  for 
"B"  to  "C"  or  "C"  to  "B"  jumps  which  would  have  a  lower  activation  energy.  Interior  "B" 
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to  "C"  to  "B"  double  jumps  are  also  possible  if  there  are  two  adjacent  vacant  sites.  It  is 
important  to  bare  in  mind  that  diffusion  effects  have  the  potential  to  dominate  the  kinetics 
of  island  growth. 

In  the  present  code,  the  crystallographic  orientations  of  surface  islands  are  defined  by  the 
substrate  coupled  with  the  matched  sizes  of  the  Ne  and  Ar  atoms.  This  gives  rise  to  a  very 
limited  number  of  boundary  orientations.  In  particular,  it  means  that  straight  boundaries 
formed  when  islands  of  different  orientation  meet  cannot  be  formed  (i.e.  boundaries  that 
are  not  simple  mirror  images  of  each  other).  However,  the  manner  in  which  such 
boundaries  develop  is  of  great  potential  importance,  because  matched  sized  systems  are 
not  the  norm.  One  way  forward  in  this  regard  is  to  use  a  coincidence  lattice  site  theoiy 
approach  where  the  boundary  is  composed  of  a  periodically  repeating  set  of  local 
environment  CA  cells.  Each  CA  cell  would  correspond  to  a  set  of  atoms  which  would 
have  evolution  or  migration  properties  which  were  predetermined  using  MD  or 
quasiharmonic  modelling  techniques.  For  visualization,  each  CA  cell  would  have  as  an 
attribute,  a  shape  which  was  a  function  of  the  nearest  neighbour  relationships  of  the 
component  atoms. 

Another  obvious  limitation  of  the  present  code  is  that  it  can  only  simulate  a  2D  surface. 
However,  it  may  not  be  too  onerous  a  task  to  extend  the  code  to  3D  if  the  aim  is  to  model 
a  deposition  process.  The  reason  is  that  the  principle  factors  in  the  kinetics  of  surface 
deposition,  and  hence  thin  layer  morphologies  are  surface  processes.  Thus  as  a  surface  is 
covered  with  new  material,  on  the  time  scale  of  surface  deposition,  atoms  now  beneath 
the  surface  may  be  regarded  as  inactive  as  far  as  the  kinetics  of  the  surface  evolution  is 
concerned.  Of  course,  they  will  continue  to  influence  the  surface  atoms  through  the 
attachment  energy.  The  obvious  way  in  which  this  could  be  achieved  in  terms  of  CA 
modelling  is  to  assign  each  CA  cell  a  height  and  a  corresponding  register  of  atoms.  Thus 
as  the  simulation  proceeds  and  the  number  of  atoms  in  the  system  increases,  the  code 
does  not  become  memory  limited  and  remains  efficient. 

Finally,  the  code  could  be  modified  to  model  the  deposition  of  molecular  species.  This 
would  be  particularly  useful  in  modelling  catalytic  processes  where  molecules  such  as  CO 
poison  surfaces.  The  essential  modification  would  be  to  associate  each  CA  cell  site  with 
an  orientation  aspect.  Also,  each  cell  may  have  a  number  of  sites  each  of  which  could 
potentially  accommodate  a  molecule,  or  part  of  a  molecule.  In  this  way,  a  variety  of 
molecular  sizes  could  be  attached  in  adjacent  cells.  This  would  have  benefits  for  the 
modelling  of  deposition  processes  where  atoms  and  molecules  in  the  incident  beam  can 
interact  with  each  other  forming  clusters  in  advance  of  deposition.  Under  such 
circumstances,  momentum  transfer  effects  may  become  important.  Although  this  is  an 
area  in  which  there  is  little  knowledge  at  present,  it  will  be  the  focus  of  a  forthcoming 
project  which  will  employ  MD  techniques  to  study  low  energy  atomic  bombardment  of 
ceramic  surfaces  (Grimes  and  Sickafiis  ICSTM  from  October  1997). 

7.3  Beyond  atomistic  based  code 

One  aspect  of  the  overall  debate  that  we  have  not  been  able  to  consider  thoroughly 
enough  is  the  extent  to  which  the  details  of  the  local  atomic  structure  influence  the  final 
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microstructure.  Certainly  we  would  expect  the  kinetics  to  be  strongly  controlled  and 
hence  the  grain  sizes  but  whether  the  type  of  microstructure  (i.e.  equiaxed ,  columnar  etc.) 
is  so  critically  affected  is  not  so  clear.  Traditionally  microstructural  evolution  models 
have  been  based  on  surface  tension  forces  and  diffusional  mass  transfer.  Both  of  these 
could  be  translated  into  CA  rule  sets.  For  example,  the  local  curvature  of  a  boundary 
which  crosses  a  CA  cell  is  a  parameter  which  would  be  effected  by  the  curvature  of 
adjacent  cells  and  would  control  mass  transport  across  the  boundary  within  the  cell. 
However,  mass  transport  is  highly  anisotropic  because  of  the  details  of  the  local  crystal 
structure.  Thus  it  is  hard  to  envisage  how  neglecting  the  effects  of  crystallography  can 
result  in  anything  other  than  an  empirically  modified  isotropic  model. 

Nevertheless,  we  cannot  remain  constrained  to  atomistic  scale  CA  cells  and  expect  to 
further  improve  computational  efficiency  to  the  extent  necessary  to  undertake  three 
dimensional  microstructural  modelling  over  long  time  scales.  Thus  it  would  seem 
inevitable  that  rather  then  associate  specific  crystallographic  positions  with  a  CA  cell,  it  is 
more  tractable  to  associate  a  crystallographic  orientation  and  thereby  introduce  anisotropy 
into  the  rule  set.  Thus,  the  probability  for  mass  transport,  evaporation  etc.  would  depend 
on  boundary  curvature  and  crystallographic  orientation.  Naturally  this  rule  set  would  have 
to  be  compiled  beforehand.  Since  experimental  data  will  usually  not  be  available,  it  will 
have  to  be  derived  using  other  modelling  techniques.  However,  the  reason  for  embarking 
on  the  present  AR-CA  study  was  partly  because  this  type  of  data  is  computationally 
expensive  to  derive  using  MD  or  energy  minimisation  techniques.  Thus  the  most  likely 
solution  would  seem  to  be  a  hierarchy  of  methods  starting  with  MD  or  energy 
minimisation,  predicting  rules  for  use  in  a  AR-CA  code  which  in  turn  provides  the 
macroscopic  rule  set. 
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8  Concluding  comments 

The  possibility  of  any  favourable  comparisons  with  microstructures  in  bulk  solids  would 
seem  to  be  unlikely  and  of  dubious  validity  at  this  stage.  Nevertheless,  the  following 
figure  is  presented  for  the  copper-phosphorus  eutectic,  an  immiscible  binary  system. 


(Micrograph  taken  from  "Materials  Science"  by  Anderson,  Leaver, 
Rawlings  &  Alexander) 


Our  present  model  assumes  no  attractive  interaction  between  "B"  type  and  "C"  type 
islands  but  only  a  very  short  range  repulsion.  In  this  respect,  the  two  types  of  gas  atom 
islands  are  acting  as  if  they  were  immiscible,  although  the  immiscibility  has  a 
crystallographic  rather  than  chemical  origin.  Furthermore,  transport  of  atoms  in  this  two 
dimensional  system  are  via  evaporation  into  the  gas  phase  rather  than  by  surface  or  solid 
state  diffusion  as  might  be  expected  in  the  three  dimensional  metallic  system. 
Nevertheless,  it  may  be  that  the  basic  model  for  a  binary  immiscibile  system  presented 
here  is  all  that  is  necessary  for  this  type  of  microstructure  to  be  formed.  The  method  of 
material  transport  may  be  largely  responsible  for  the  kinetics  of  growth  rather  than  for  the 
gross  microstructural  features.  Therefore,  although  the  authors  wish  to  stress  that 
considerably  more  work  is  necessary  before  firm  conclusions  can  be  drawn,  (in  particular 
diffusion  rules  must  be  introduced  into  the  model)  these  results  are  certainly  encouraging. 
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